## wolves and sheeps model example 2
# same values from NetLogo example

# notice D is diffusion rate. eg: rate/(l^2)

volterra.setup <- function(W, H, X0, C, D)
{
	# R1: X1    -> 2.X1
	# R2: X1+X2 -> 2.X2
	# R3: X2    -> 0
	transitions <- function(X) c(
		C[1]*X[1],
		C[2]*X[1]*X[2],
		C[3]*X[2])
	V <- matrix(c(  # stoichiometry matrix
		+1, 0,
		-1,+1,
		 0,-1), ncol=2, byrow=TRUE)

	X <- matrix(rep(0,W*H*2), nrow=W*H)
	for(j in 1:2)
		for(i in trunc(runif(X0[j],1,W*H+1)))
			X[i,j] <- X[i,j]+1

	source('spatial.R')
	ctx <- new.env()
	spatial.setup(ctx, X, transitions, V, W, H, D)
	ctx$ts <- NULL
	ctx$xs <- NULL
	ctx$ys <- NULL
	ctx
}

volterra.step <- function(ctx)
{
	oldt <- ctx$t
	ctx$ts <- c(ctx$ts, ctx$t)
	ctx$xs <- c(ctx$xs, sum(ctx$X[,1]))
	ctx$ys <- c(ctx$ys, sum(ctx$X[,2]))
	spatial.step(ctx)
	if(ctx$t < oldt)  # HACK
		ctx$t <- oldt
}

volterra.plot <- function(ctx)
{
	dev.hold()
	with(ctx, {
		Y <- ifelse(X[,1]>0,1,0) + ifelse(X[,2]>0,2,0)
		par(mar=c(0,0,0,0))
		image(matrix(Y,W,H), col=c('green','white','black','red'), zlim=c(0,3), axes=FALSE)
		par(mar=c(2,2,0,0))
		plot(ts, xs, type='l', xlim=c(0,xmax), ylim=c(0,max(xs,ys,1000)), lwd=2, col='blue')
		lines(ts, ys, lwd=2, col='red')
	})
	dev.flush()
}

volterra.run <- function(W, H, X0, C, D, tmax=0)
{
	ctx <- volterra.setup(W, H, X0, C, D)
	ctx$xmax <- ifelse(tmax==0,50,tmax)

	if(tmax == 0) {
		tryCatch(dev.off(), error=function(e) NULL)  # close window
		dev.new(width=10, height=5)
	}
	par(mfrow=c(1,2))

	i <- 0
	while(ctx$t < ctx$xmax) {
		volterra.step(ctx)
		if(tmax == 0) {
			if(i %% 50 == 0)
				volterra.plot(ctx)
			i <- i+1
		}
	}
	if(tmax > 0)
		volterra.plot(ctx)
}

#volterra.run(50, 50, c(100, 50), c(0.2,0.5,0.1), c(1,5)/9)
#volterra.run(1, 1, c(100, 50), c(1,0.005,0.6), c(0,0))
